module cubed_sphere_mesh_mod

  use mpi
  use kinds_mod
  use mesh_const_mod
  use cubed_sphere_math_mod
  use cubed_sphere_domain_mod

  implicit none

  private

  public cubed_sphere_domain_type
  public cubed_sphere_mesh_type
  public equiangular

  type cubed_sphere_mesh_type
    logical  :: initialized = .false.
    integer  :: pid         = MPI_PROC_NULL
    integer  :: proj_type   = equiangular
    integer  :: domain_type = global_domain
    integer  :: ndom        = 6      ! Number of domains
    integer  :: dms         = 0      ! Domain start index
    integer  :: dme         = 0      ! Domain end index
    real(r8) :: rlon        = 0      ! Longitude of the rotated origin point
    real(r8) :: rlat        = 0      ! Latitude of the rotated origin point
    type(cubed_sphere_domain_type), allocatable :: domain(:)
    ! Procedures
    procedure(locate_domain_interface           ), pointer, nopass :: locate_domain
    procedure(sphere_to_cube_interface          ), pointer, nopass :: sphere_to_cube
    procedure(cube_to_sphere_interface          ), pointer, nopass :: cube_to_sphere
    procedure(translate_idx_along_edge_interface), pointer,   pass :: translate_idx_along_edge => null()
  contains
    procedure :: init             => cubed_sphere_mesh_init
    procedure :: set_hor_metrics  => cubed_sphere_mesh_set_horizontal_metrics
    procedure :: set_ghost_points => cubed_sphere_mesh_set_ghost_points
    procedure :: get_dom_Jh       => cubed_sphere_mesh_get_dom_Jh
    procedure :: get_dom_Gh       => cubed_sphere_mesh_get_dom_Gh
    procedure :: get_dom_iGh      => cubed_sphere_mesh_get_dom_iGh
    procedure :: get_dom_iG       => cubed_sphere_mesh_get_dom_iG
    procedure :: get_dom_Ah       => cubed_sphere_mesh_get_dom_Ah
    procedure :: get_dom_iAh      => cubed_sphere_mesh_get_dom_iAh
    procedure :: get_dom_CS       => cubed_sphere_mesh_get_dom_CS
    procedure :: get_gst_Jh       => cubed_sphere_mesh_get_gst_Jh
    procedure :: get_gst_Gh       => cubed_sphere_mesh_get_gst_Gh
    procedure :: get_gst_iGh      => cubed_sphere_mesh_get_gst_iGh
    procedure :: get_gst_Ah       => cubed_sphere_mesh_get_gst_Ah
    procedure :: get_gst_iAh      => cubed_sphere_mesh_get_gst_iAh
    procedure :: clear            => cubed_sphere_mesh_clear
    procedure :: get_params       => cubed_sphere_mesh_get_params
    final :: cubed_sphere_mesh_final
  end type cubed_sphere_mesh_type
  
  interface
    subroutine translate_idx_along_edge_interface(this, idom, orient, ic, jc, ipe, ic_ngb, jc_ngb, ipe_ngb)
      import cubed_sphere_mesh_type
      class(cubed_sphere_mesh_type), intent(in) :: this
      integer, intent(in ) :: idom    ! Domain index
      integer, intent(in ) :: orient  ! Edge orient
      integer, intent(in ) :: ic      ! Cell index along x-axis
      integer, intent(in ) :: jc      ! Cell index along y-axis
      integer, intent(in ) :: ipe     ! Edge point index
      integer, intent(out) :: ic_ngb  ! Adjacent cell index along x-axis
      integer, intent(out) :: jc_ngb  ! Adjacent cell index along y-axis
      integer, intent(out) :: ipe_ngb ! Adjacent edge point index
    
    end subroutine translate_idx_along_edge_interface
  end interface

contains

  subroutine cubed_sphere_mesh_init(this, proj_type, hdim, vdim, nc, ndx, ndy, idom, hw, rw, rwv, neq, xeq, weq, &
                                    ids, jds, nbx, nby, nz, pid, cart_dims, cart_coords, stage)

    class(cubed_sphere_mesh_type), intent(inout) :: this
    integer , intent(in)           :: proj_type      ! Projection type (only equiangular for now)
    integer , intent(in), optional :: hdim           ! horizontal reconstruction dimension, choose from 2 or 3
    integer , intent(in), optional :: vdim           ! vertical reconstruction dimension, choose from 1 to 3
    integer , intent(in), optional :: nc             ! Resolution parameter
    integer , intent(in), optional :: ndx            ! Number of domain cells along x-axis
    integer , intent(in), optional :: ndy            ! Number of domain cells along y-axis
    integer , intent(in), optional :: idom           ! Index for active domain
    integer , intent(in), optional :: hw             ! Halo width
    integer , intent(in), optional :: rw             ! Horizontal reconstruction width
    integer , intent(in), optional :: rwv            ! Vertical reconstruction width
    integer , intent(in), optional :: neq            ! Number of quadrature points along cell edge
    real(16), intent(in), optional :: xeq(:)         ! Positions of quadrature points along cell edge
    real(r8), intent(in), optional :: weq(:)         ! Weights of quadrature points along cell edge
    integer , intent(in), optional :: ids            ! Start index along x axis
    integer , intent(in), optional :: jds            ! Start index along y axis
    integer , intent(in), optional :: nbx            ! Number of grids along x axis
    integer , intent(in), optional :: nby            ! Number of grids along y axis
    integer , intent(in), optional :: nz             ! Number of vertical levels
    integer , intent(in), optional :: pid            ! Process ID
    integer , intent(in), optional :: cart_dims(2)   ! Process Cartesian decomposition dimensions
    integer , intent(in), optional :: cart_coords(2) ! Process Cartesian decomposition coordinates
    integer , intent(in), optional :: stage          ! Initialization stage (there are 2 stages)

    integer i

    if (merge(stage, 1, present(stage)) == 1) then
      call this%clear()
      allocate(this%domain(6))
    end if

    if (present(pid)) this%pid = pid
    this%proj_type = proj_type

    ! Domain topology
    !        _____
    !       |     |
    !       |  5  |
    !  _____|_____|___________
    ! |     |     |     |     |
    ! |  4  |  1  |  2  |  3  |
    ! |_____|_____|_____|_____|
    !       |     |
    !       |  6  |
    !       |_____|
    !
    call this%domain(1)%connect(1, this%domain(4), this%domain(2), this%domain(5), this%domain(6))
    call this%domain(2)%connect(2, this%domain(1), this%domain(3), this%domain(5), this%domain(6))
    call this%domain(3)%connect(3, this%domain(2), this%domain(4), this%domain(5), this%domain(6))
    call this%domain(4)%connect(4, this%domain(3), this%domain(1), this%domain(5), this%domain(6))
    call this%domain(5)%connect(5, this%domain(4), this%domain(2), this%domain(3), this%domain(1))
    call this%domain(6)%connect(6, this%domain(4), this%domain(2), this%domain(1), this%domain(3))

    ! Set neighbor domain if necessary.
    if (present(idom)) then
      call this%domain(idom)%init(proj_type, hdim=hdim, vdim=vdim, nc=nc, ndx=ndx, ndy=ndy, hw=hw, rw=rw, rwv=rwv, &
                                  neq=neq, xeq=xeq, weq=weq, ids=ids, jds=jds, nbx=nbx, nby=nby, &
                                  nz=nz, active=.true., pid=pid, cart_dims=cart_dims, cart_coords=cart_coords)
      this%domain(idom)%active = .true.
      this%dms = idom
      this%dme = idom
    else
      do i = 1, this%ndom
        call this%domain(i)%init(proj_type, hdim=hdim, vdim=vdim, nc=nc, ndx=ndx, ndy=ndy, hw=hw, rw=rw, rwv=rwv, &
                                 neq=neq, xeq=xeq, weq=weq, nz=nz, pid=pid, cart_dims=cart_dims, cart_coords=cart_coords)
        this%domain(i)%active = .true.
      end do
      this%dms = 1
      this%dme = this%ndom
    end if

    select case (proj_type)
    case (equiangular)
      this%locate_domain  => equiangular_locate_domain
      this%sphere_to_cube => equiangular_sphere_to_cube
      this%cube_to_sphere => equiangular_cube_to_sphere
    end select
    
    if(hdim==2)then
      this%translate_idx_along_edge => translate_idx_along_edge_2d
    elseif(hdim==3)then
      this%translate_idx_along_edge => translate_idx_along_edge_3d
    endif

    this%initialized = .true.

  end subroutine cubed_sphere_mesh_init

  subroutine cubed_sphere_mesh_set_horizontal_metrics(this)

    class(cubed_sphere_mesh_type), intent(inout) :: this

    integer idom

    do idom = 1, this%ndom
      if (this%domain(idom)%initialized) call this%domain(idom)%set_hor_metrics()
    end do

  end subroutine cubed_sphere_mesh_set_horizontal_metrics

  subroutine cubed_sphere_mesh_set_ghost_points(this)

    class(cubed_sphere_mesh_type), intent(inout) :: this

    integer idom

    do idom = 1, this%ndom
      if (this%domain(idom)%initialized) call this%domain(idom)%count_ghost_points()
    end do

    do idom = 1, this%ndom
      if (this%domain(idom)%initialized) call this%domain(idom)%allocate_ghost_points()
    end do

    do idom = 1, this%ndom
      if (this%domain(idom)%initialized) call this%domain(idom)%add_ghost_points()
    end do

    do idom = 1, this%ndom
      if (this%domain(idom)%initialized) call this%domain(idom)%set_ghost_ij()
    end do

    do idom = 1, this%ndom
      if (this%domain(idom)%active) then
        this%domain(idom)%ngb(1)%ex_depth = maxval(this%domain(idom)%ngb(1)%ptr%ngb(:)%ex_depth)
        this%domain(idom)%ngb(2)%ex_depth = maxval(this%domain(idom)%ngb(2)%ptr%ngb(:)%ex_depth)
        this%domain(idom)%ngb(3)%ex_depth = maxval(this%domain(idom)%ngb(3)%ptr%ngb(:)%ex_depth)
        this%domain(idom)%ngb(4)%ex_depth = maxval(this%domain(idom)%ngb(4)%ptr%ngb(:)%ex_depth)
      end if
    end do

  end subroutine cubed_sphere_mesh_set_ghost_points

  subroutine cubed_sphere_mesh_clear(this)

    class(cubed_sphere_mesh_type), intent(inout) :: this

    if (allocated(this%domain)) deallocate(this%domain)

    this%initialized = .false.

  end subroutine cubed_sphere_mesh_clear

  subroutine translate_idx_along_edge_2d(this, idom, orient, ic, jc, ipe, ic_ngb, jc_ngb, ipe_ngb)

    class(cubed_sphere_mesh_type), intent(in) :: this
    integer, intent(in ) :: idom    ! Domain index
    integer, intent(in ) :: orient  ! Edge orient
    integer, intent(in ) :: ic      ! Cell index along x-axis
    integer, intent(in ) :: jc      ! Cell index along y-axis
    integer, intent(in ) :: ipe     ! Edge point index
    integer, intent(out) :: ic_ngb  ! Adjacent cell index along x-axis
    integer, intent(out) :: jc_ngb  ! Adjacent cell index along y-axis
    integer, intent(out) :: ipe_ngb ! Adjacent edge point index

    integer i, ipe1, neq, nc, pve
    
    pve  = this%domain(idom)%pve
    ipe1 = ipe - pve
    neq  = this%domain(idom)%neq
    nc   = this%domain(idom)%nc

    select case (orient)
    case (left)
      i = ipe1 - 3 * neq
      select case (this%domain(idom)%ngb(1)%code)
      case (14, 21, 32, 43)
        ic_ngb  = nc
        jc_ngb  = jc
        ipe_ngb = i + neq
      case (64)
        ic_ngb  = jc
        jc_ngb  = 1
        ipe_ngb = i
      case (54)
        ic_ngb  = nc - jc + 1
        jc_ngb  = nc
        ipe_ngb = 3 * neq - i + 1
      end select
    case (right)
      i = ipe1 - 1 * neq
      select case (this%domain(idom)%ngb(2)%code)
      case (12, 23, 34, 41)
        ic_ngb  = 1
        jc_ngb  = jc
        ipe_ngb = i + 3 * neq
      case (52)
        ic_ngb  = jc
        jc_ngb  = nc
        ipe_ngb = i + 2 * neq
      case (62)
        ic_ngb  = nc - jc + 1
        jc_ngb  = 1
        ipe_ngb = neq - i + 1
      end select
    case (bottom)
      i = ipe1
      select case (this%domain(idom)%ngb(3)%code)
      case (16, 51)
        ic_ngb  = ic
        jc_ngb  = nc
        ipe_ngb = i + 2 * neq
      case (26)
        ic_ngb  = nc
        jc_ngb  = nc - ic + 1
        ipe_ngb = 2 * neq - i + 1
      case (46)
        ic_ngb  = 1
        jc_ngb  = ic
        ipe_ngb = i + 3 * neq
      case (36, 63)
        ic_ngb  = nc - ic + 1
        jc_ngb  = 1
        ipe_ngb = neq - i + 1
      end select
    case (top)
      i = ipe1 - 2 * neq
      select case (this%domain(idom)%ngb(4)%code)
      case (15, 61)
        ic_ngb  = ic
        jc_ngb  = 1
        ipe_ngb = i
      case (45)
        ic_ngb  = 1
        jc_ngb  = nc - ic + 1
        ipe_ngb = 4 * neq - i + 1
      case (25)
        ic_ngb  = nc
        jc_ngb  = ic
        ipe_ngb = i + neq
      case (35, 53)
        ic_ngb  = nc - ic + 1
        jc_ngb  = nc
        ipe_ngb = 3 * neq - i + 1
      end select
    end select

    ipe_ngb = ipe_ngb + pve

  end subroutine translate_idx_along_edge_2d

  subroutine translate_idx_along_edge_3d(this, idom, orient, ic, jc, ipe, ic_ngb, jc_ngb, ipe_ngb)

    class(cubed_sphere_mesh_type), intent(in) :: this
    integer, intent(in ) :: idom    ! Domain index
    integer, intent(in ) :: orient  ! Edge orient
    integer, intent(in ) :: ic      ! Cell index along x-axis
    integer, intent(in ) :: jc      ! Cell index along y-axis
    integer, intent(in ) :: ipe     ! Edge point index
    integer, intent(out) :: ic_ngb  ! Adjacent cell index along x-axis
    integer, intent(out) :: jc_ngb  ! Adjacent cell index along y-axis
    integer, intent(out) :: ipe_ngb ! Adjacent edge point index

    integer i, k, ipe1, neq, nfq, nc, pve
    
    pve  = this%domain(idom)%pve
    ipe1 = ipe - pve
    neq  = this%domain(idom)%neq
    nfq  = this%domain(idom)%nfq
    nc   = this%domain(idom)%nc

    select case (orient)
    case (left)
      i = ipe1 - 3 * nfq
      select case (this%domain(idom)%ngb(1)%code)
      case (14, 21, 32, 43)
        ic_ngb  = nc
        jc_ngb  = jc
        ipe_ngb = i + nfq
      case (64)
        ic_ngb  = jc
        jc_ngb  = 1
        ipe_ngb = i
      case (54)
        ic_ngb  = nc - jc + 1
        jc_ngb  = nc
        k = ( i - 1 ) / neq + 1
        ipe_ngb = 2 * nfq + k * neq - mod( i - 1, neq )
      end select
    case (right)
      i = ipe1 - 1 * nfq
      select case (this%domain(idom)%ngb(2)%code)
      case (12, 23, 34, 41)
        ic_ngb  = 1
        jc_ngb  = jc
        ipe_ngb = i + 3 * nfq
      case (52)
        ic_ngb  = jc
        jc_ngb  = nc
        ipe_ngb = i + 2 * nfq
      case (62)
        ic_ngb  = nc - jc + 1
        jc_ngb  = 1
        k = ( i - 1 ) / neq + 1
        ipe_ngb = k * neq - mod( i - 1, neq )
      end select
    case (bottom)
      i = ipe1
      select case (this%domain(idom)%ngb(3)%code)
      case (16, 51)
        ic_ngb  = ic
        jc_ngb  = nc
        ipe_ngb = i + 2 * nfq
      case (26)
        ic_ngb  = nc
        jc_ngb  = nc - ic + 1
        k = ( i - 1 ) / neq + 1
        ipe_ngb = nfq + k * neq - mod( i - 1, neq )
      case (46)
        ic_ngb  = 1
        jc_ngb  = ic
        ipe_ngb = i + 3 * nfq
      case (36, 63)
        ic_ngb  = nc - ic + 1
        jc_ngb  = 1
        k = ( i - 1 ) / neq + 1
        ipe_ngb = k * neq - mod( i - 1, neq )
      end select
    case (top)
      i = ipe1 - 2 * nfq
      select case (this%domain(idom)%ngb(4)%code)
      case (15, 61)
        ic_ngb  = ic
        jc_ngb  = 1
        ipe_ngb = i
      case (45)
        ic_ngb  = 1
        jc_ngb  = nc - ic + 1
        k = ( i - 1 ) / neq + 1
        ipe_ngb = 3 * nfq + k * neq - mod( i - 1, neq )
      case (25)
        ic_ngb  = nc
        jc_ngb  = ic
        ipe_ngb = i + nfq
      case (35, 53)
        ic_ngb  = nc - ic + 1
        jc_ngb  = nc
        k = ( i - 1 ) / neq + 1
        ipe_ngb = 2 * nfq + k * neq - mod( i - 1, neq )
      end select
    end select

    ipe_ngb = ipe_ngb + pve

  end subroutine translate_idx_along_edge_3d
  
  function cubed_sphere_mesh_get_dom_Jh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res
    integer :: idom,p,i,j,k
    
    associate(Jh => this%domain(idom)%Jh, &
              r  => this%domain(idom)%r   )
      res = Jh(p,i,j) * r(p,i,j,k) * r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_Jh
  
  function cubed_sphere_mesh_get_dom_Ah(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(Ah => this%domain(idom)%Ah, &
              r  => this%domain(idom)%r   )
      res = Ah(:,:,p,i,j) * r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_Ah
  
  function cubed_sphere_mesh_get_dom_iAh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(iAh => this%domain(idom)%iAh, &
              r   => this%domain(idom)%r    )
      res = iAh(:,:,p,i,j) / r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_iAh
  
  function cubed_sphere_mesh_get_dom_CS(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(3,3,3)
    integer :: idom,p,i,j,k
    
    associate(CS => this%domain(idom)%CS, &
              r   => this%domain(idom)%r  )
      res        = CS(:,:,:,p,i,j)
      res(1,3,1) = CS(1,3,1,p,i,j) / r(p,i,j,k)
      res(3,1,1) = CS(3,1,1,p,i,j) / r(p,i,j,k)
      res(2,3,2) = CS(2,3,2,p,i,j) / r(p,i,j,k)
      res(3,2,2) = CS(3,2,2,p,i,j) / r(p,i,j,k)
      res(1,1,3) = CS(1,1,3,p,i,j) * r(p,i,j,k)
      res(1,2,3) = CS(1,2,3,p,i,j) * r(p,i,j,k)
      res(2,1,3) = CS(2,1,3,p,i,j) * r(p,i,j,k)
      res(2,2,3) = CS(2,2,3,p,i,j) * r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_CS
  
  function cubed_sphere_mesh_get_dom_Gh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(Gh => this%domain(idom)%Gh, &
              r  => this%domain(idom)%r   )
      res = Gh(:,:,p,i,j) * r(p,i,j,k) * r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_Gh
  
  function cubed_sphere_mesh_get_dom_iGh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(iGh => this%domain(idom)%iGh, &
              r   => this%domain(idom)%r    )
      res = iGh(:,:,p,i,j) / r(p,i,j,k) / r(p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_iGh
  
  function cubed_sphere_mesh_get_dom_iG(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(3,3)
    integer  :: idom,p,i,j,k
    real(r8) :: iGh(2,2)
    
    iGh = this%get_dom_iGh(idom,p,i,j,k)
    
    associate(Jv  => this%domain(idom)%Jv,  &
              iGv => this%domain(idom)%iGv, &
              iG  => this%domain(idom)%iG   )
      ! res          = 0
      ! res(1:2,1:2) = iGh(1:2,1:2)
      ! res(1  ,3  ) = iGv(1,3,p,i,j,k) * iGh(1,1) + iGv(2,3,p,i,j,k) * iGh(1,2)
      ! res(2  ,3  ) = iGv(1,3,p,i,j,k) * iGh(2,1) + iGv(2,3,p,i,j,k) * iGh(2,2)
      ! res(3  ,3  ) = 1._r_kind / Jv(p,i,j,k)

      ! res          = 0
      ! res(1:2,1:2) = iGh(1:2,1:2)
      ! res(1  ,3  ) = iG(1,3,p,i,j,k)
      ! res(2  ,3  ) = iG(2,3,p,i,j,k)
      ! res(3  ,3  ) = 1._r_kind / Jv(p,i,j,k)

      res          = 0
      res(1:2,1:2) = iGh(1:2,1:2)
      res(1  ,3  ) = iG(1,3,p,i,j,k)
      res(2  ,3  ) = iG(2,3,p,i,j,k)
      res(3  ,3  ) = iG(3,3,p,i,j,k)
    end associate
  end function cubed_sphere_mesh_get_dom_iG
  
  function cubed_sphere_mesh_get_gst_Jh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res
    integer :: idom,p,i,j,k
    
    associate(Jh => this%domain(idom)%ghost(i,j)%Jh, &
              r  => this%domain(idom)%ghost(i,j)%r   )
      res = Jh(p) * r(p,k) * r(p,k)
    end associate
  end function cubed_sphere_mesh_get_gst_Jh
  
  function cubed_sphere_mesh_get_gst_Ah(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(Ah => this%domain(idom)%ghost(i,j)%Ah, &
              r  => this%domain(idom)%ghost(i,j)%r   )
      res = Ah(:,:,p) * r(p,k)
    end associate
  end function cubed_sphere_mesh_get_gst_Ah
  
  function cubed_sphere_mesh_get_gst_iAh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(iAh => this%domain(idom)%ghost(i,j)%iAh, &
              r   => this%domain(idom)%ghost(i,j)%r    )
      res = iAh(:,:,p) / r(p,k)
    end associate
  end function cubed_sphere_mesh_get_gst_iAh
  
  function cubed_sphere_mesh_get_gst_Gh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(Gh => this%domain(idom)%ghost(i,j)%Gh, &
              r  => this%domain(idom)%ghost(i,j)%r   )
      res = Gh(:,:,p) * r(p,k) * r(p,k)
    end associate
  end function cubed_sphere_mesh_get_gst_Gh
  
  function cubed_sphere_mesh_get_gst_iGh(this,idom,p,i,j,k) result (res)
    class(cubed_sphere_mesh_type), intent(inout) :: this
    real(r8) :: res(2,2)
    integer :: idom,p,i,j,k
    
    associate(iGh => this%domain(idom)%ghost(i,j)%iGh, &
              r   => this%domain(idom)%ghost(i,j)%r    )
      res = iGh(:,:,p) / r(p,k) / r(p,k)
    end associate
  end function cubed_sphere_mesh_get_gst_iGh
  
  subroutine cubed_sphere_mesh_get_params(this, idom,             &
                                          nx, ny, nz,             &
                                          nbx, nby, ndx, ndy,     &
                                          ims, ime, ids, ide,     &
                                          jms, jme, jds, jde,     &
                                          kms, kme, kds, kde,     &
                                          dms, dme,               &
                                          ipds, ipde, jpds, jpde, &
                                          irs, ire, jrs, jre,     &
                                          its, ite, jts, jte,     &
                                          igs, ige, jgs, jge,     &
                                          npt, nv, ncq, neq,      &
                                          nfq, nvq, nfqv,         &
                                          pc, pvs, pve, pes, pee, &
                                          pqs, pqe, pgs, pge,     &
                                          xeq, weq1d, weq2d,      &
                                          max_ngp)

    class(cubed_sphere_mesh_type), intent(in) :: this
    integer , intent(in ), optional :: idom
    integer , intent(out), optional :: nx, ny, nz
    integer , intent(out), optional :: nbx, nby
    integer , intent(out), optional :: ndx, ndy
    integer , intent(out), optional :: ims, ime, ids, ide
    integer , intent(out), optional :: jms, jme, jds, jde
    integer , intent(out), optional :: kms, kme, kds, kde
    integer , intent(out), optional :: dms, dme
    integer , intent(out), optional :: ipds, ipde, jpds, jpde
    integer , intent(out), optional :: irs, ire, jrs, jre
    integer , intent(out), optional :: its, ite, jts, jte
    integer , intent(out), optional :: igs, ige, jgs, jge
    integer , intent(out), optional :: npt
    integer , intent(out), optional :: nv
    integer , intent(out), optional :: ncq
    integer , intent(out), optional :: neq
    integer , intent(out), optional :: nfq
    integer , intent(out), optional :: nvq
    integer , intent(out), optional :: nfqv
    integer , intent(out), optional :: pc
    integer , intent(out), optional :: pvs
    integer , intent(out), optional :: pve
    integer , intent(out), optional :: pes(6)
    integer , intent(out), optional :: pee(6)
    integer , intent(out), optional :: pqs
    integer , intent(out), optional :: pqe
    integer , intent(out), optional :: pgs
    integer , intent(out), optional :: pge
    real(16), intent(out), optional :: xeq(:)
    real(r8), intent(out), optional :: weq1d(:)
    real(r8), intent(out), optional :: weq2d(:)
    integer , intent(out), optional :: max_ngp
    integer i

    if (present(idom)) then
      i = idom
    else
      i = this%dms
    end if

    if (present(nx )) nx  = this%domain(i)%nx
    if (present(ny )) ny  = this%domain(i)%ny
    if (present(nz )) nz  = this%domain(i)%nz
    if (present(nbx)) nbx = this%domain(i)%nbx
    if (present(nby)) nby = this%domain(i)%nby
    if (present(ndx)) ndx = this%domain(i)%ndx
    if (present(ndy)) ndy = this%domain(i)%ndy
    if (present(ims)) ims = this%domain(i)%ims
    if (present(ime)) ime = this%domain(i)%ime
    if (present(ids)) ids = this%domain(i)%ids
    if (present(ide)) ide = this%domain(i)%ide
    if (present(jms)) jms = this%domain(i)%jms
    if (present(jme)) jme = this%domain(i)%jme
    if (present(jds)) jds = this%domain(i)%jds
    if (present(jde)) jde = this%domain(i)%jde
    if (present(kms)) kms = this%domain(i)%kms
    if (present(kme)) kme = this%domain(i)%kme
    if (present(kds)) kds = this%domain(i)%kds
    if (present(kde)) kde = this%domain(i)%kde
    if (present(dms)) dms = this%dms
    if (present(dme)) dme = this%dme

    if (present(ipds)) ipds = this%domain(i)%ipds
    if (present(ipde)) ipde = this%domain(i)%ipde
    if (present(jpds)) jpds = this%domain(i)%jpds
    if (present(jpde)) jpde = this%domain(i)%jpde
    
    if (present(irs)) irs = this%domain(i)%irs
    if (present(ire)) ire = this%domain(i)%ire
    if (present(jrs)) jrs = this%domain(i)%jrs
    if (present(jre)) jre = this%domain(i)%jre
    
    if (present(its)) its = this%domain(i)%its
    if (present(ite)) ite = this%domain(i)%ite
    if (present(jts)) jts = this%domain(i)%jts
    if (present(jte)) jte = this%domain(i)%jte
    
    if (present(igs)) igs = this%domain(i)%igs
    if (present(ige)) ige = this%domain(i)%ige
    if (present(jgs)) jgs = this%domain(i)%jgs
    if (present(jge)) jge = this%domain(i)%jge

    if (present(npt )) npt  = this%domain(i)%npt
    if (present(nv  )) nv   = this%domain(i)%nv
    if (present(ncq )) ncq  = this%domain(i)%ncq
    if (present(neq )) neq  = this%domain(i)%neq
    if (present(nfqv)) nfqv = this%domain(i)%nfqv
    if (present(nfq )) nfq  = this%domain(i)%nfq
    if (present(nvq )) nvq  = this%domain(i)%nvq
    if (present(pc  )) pc   = this%domain(i)%pc
    if (present(pvs )) pvs  = this%domain(i)%pvs
    if (present(pve )) pve  = this%domain(i)%pve
    if (present(pes )) pes  = this%domain(i)%pes
    if (present(pee )) pee  = this%domain(i)%pee
    if (present(pqs )) pqs  = this%domain(i)%pqs
    if (present(pqe )) pqe  = this%domain(i)%pqe
    if (present(pgs )) pgs  = this%domain(i)%pgs
    if (present(pge )) pge  = this%domain(i)%pge

    if (present(xeq  )) xeq   = this%domain(i)%xeq
    if (present(weq1d)) weq1d = this%domain(i)%weq1d
    if (present(weq2d)) weq2d = this%domain(i)%weq2d

    if (present(max_ngp)) max_ngp = this%domain(i)%max_ngp

  end subroutine cubed_sphere_mesh_get_params

  subroutine cubed_sphere_mesh_final(this)

    type(cubed_sphere_mesh_type), intent(inout) :: this

    call this%clear()

  end subroutine cubed_sphere_mesh_final

end module cubed_sphere_mesh_mod
